Решить задачу квадратичного программирования $$ \begin{align}{c} \min -\bar{p}^Tx + \mu x^T \Sigma x\\ \sum x_i = 1,\ x_i \geq 0 \end{align} $$ для большого количества различных положительных $\mu$ (например, 100 значений, логарифмически распределенных от 1 до $10^7$). Нарисовать графики $\bar{p}^T x$ от стандартного отклонения $(x^T \sigma x)^{1/2}$, а также картинки, как на (рис. 4.12, стр. 187, Boyd).
В наборе четыре инструмента, один из которых безрисковый и дает доход в 3%.
Для использования этого метода задача переформулируется следующим образом: $$ \begin{array}{rl} \text{minimize} & (m / \varepsilon) f_0(x) + \phi(x)\\ \text{subject to} & Ax = b \end{array} $$
Идея метода состоит в решении последовательности минимизационных подзадач: вычисляется значение $x^*(t)$ для возрастающей последовательности из $t$, пока $t < m / \varepsilon$.
Алгоритм выглядит следующим образом: сначала задаются $x, t := t^{(0)} = 0, \mu > 1, \varepsilon > 0$. Затем повторяются следующие операции:
Метод останавливается, когда $f_0(x) - p^* \leq \varepsilon$ (следует из $f_0(x^(t)) - p^ \leq m / t))
Центрирование обычно делается методом Ньютона, начиная от текущего значения $x$.
Выбор $\mu$ подразумевает компромисс: большое значение $\mu$ означает меньшее число внешних итераций, большее -- большее; обычно выбирают значения от 10 до 20.
Используется несколько эвристик для выбора $t^{(0)}$.
Число внешних итераций (центрирования) задается, как $$ \bigg\lceil \dfrac{\log(m / \varepsilon t^{(0)}))}{\log \mu}\bigg\rceil $$ плюс изначальный шаг центрирования (для вычисления $x^*(t^{(0)})$).
Задача центрирования: $$ \min t f_0(x) + \phi(x) $$
Для этого необходимо иметь результаты анализа сходимости метода Ньютона.
Итеративный метод поиска корней уравнения $f'(x) = 0$, при этом $f$ должна быть дифференцируемой. $$ x_{n+1} = x_n - \dfrac{f'(x_n)}{f''(x_n)} $$
|
![]() |
|
![]() |
Предполагается, что дана нормально распределенная случайная величина $p\sim\mathcal{N}(\bar{p}, \Sigma)$. Требуется сформулировать задачу $$ \begin{array}{c} \max \bar{p}^T x\\ prob(p^Tx\leq 0) \leq \eta\\ \sum x_i = 1,\ x_i \geq 0 \end{array} $$ как задачу выпуклой оптимизации с параметром $\eta < 1/2$, и решить ее для большого числа разных $\eta$ от $10^{-4}$ до $10^{-1}$, нарисовав оптимальные $\bar{p}^T x$ от $\eta$.
Но $$ prob(p^Tx\leq 0) = \Phi\left( \dfrac{-\bar{p}^T x}{\lVert \Sigma^{1/2} x \rVert_2} \right) = 1 - \Phi\left( \dfrac{\bar{p}^T x}{\lVert \Sigma^{1/2} x \rVert_2} \right) $$ Теперь $$ \begin{array}{c} \max \bar{p}^T x\\ -\bar{p}^Tx + \Phi^{-1}(1 - \eta)\lVert \Sigma^{1/2} x \rVert_2 \leq 0\\ \sum x_i = 1,\ x_i \geq 0 \end{array} $$
То, что мы вывели раньше, формулируется в терминах SOCP: $$ \begin{array}{rl} \text{minimize} & f^T x\\ \text{subject to} & -\lVert A_i x + b_i \rVert_2 \leq c_i^T x + d_i, i = 1,\ldots,m\\ \end{array} $$ Сформулируем нашу задачу в этих терминах $$ \begin{array}{rl} \text{minimize} & -\bar{p}^T x\\ \text{subject to} & \lVert \Sigma^{1/2} x \rVert_2 \leq \dfrac{1}{\Phi^{-1}(1 - \eta)}\bar{p}^Tx \\ & \sum x_i = 1\\ & \mathbf{diag}(x_i) \succeq 0 \end{array} $$ Ограничения встроились бы в виде стандартной комбинации двух неравенств, которую мы бы преобразовали так, чтобы все выражалось через норму, но здесь этого делать не будем, так как мы найдем новую эквивалентную задачу.
Задача такого вида, где $F_i, G \in \mathbb{S}^k$, называется задачей полуопределенного программирования, и именно ее мы и будем решать численно. Для задачи конического программирования второго порядка $$ \begin{array}{rl} \text{minimize} & f^T x\\ \text{subject to} & \lVert A_i x + b_i \rVert_2 \leq c_i^T x + d_i, i = 1,\ldots,m\\ \end{array} $$ существует эквивалентная задача полуопределенного программирования: $$ \begin{array}{rl} \text{minimize} & f^T x\\ \text{subject to} & \left[ \begin{array}{cc} (c_i^T x + d_i)I & A_ix+b_i\\ (A_i x + b_i)^T & c_i^T x + d_i \end{array} \right] \succeq 0,\ i = 1,\ldots,m \end{array} $$ В нее уже можно будет внести наши ограничения явно.
Ограничение на единичную сумму иксов введем, разложив покомпонентно по иксам приведенную выше матрицу, а затем приписав к каждой из матриц справа снизу число 1, и определив соответствующий элемент матрицы $G$, как 1.
Как мы будем разделять эту матрицу на компоненты? Подставляя в нее в качестве $x$ элементы базиса (векторы с единственной ненулевой координатой, равной единице).
Мы использовали метод CSDP.
library(Rcsdp)
library(expm)
library(plotly)
p_bar <- c(0.12, 0.10, 0.07, 0.03)
Sigma <- matrix(
c(
0.0064, 0.0008, -0.0011, 0,
0.0008, 0.0025, 0, 0,
-0.0011, 0, 0.0004, 0,
0, 0, 0, 0
),
nrow = 4
)
Sigma_sqrt <- sqrtm(Sigma)
optimal_res <- sapply(
seq(0.0001, 0.1, 0.0001),↔
function(eta) {
p_bar.mat <- matrix(p_bar)
C <- list(
matrix(0, length(p_bar) + 1, length(p_bar)+ 1),
c(rep(0, length(p_bar)), 1, -1)
)
A <- lapply(↔
)
b <- (-1) * p_bar
K <- list(type = c("s", "l"), size = c(length(p_bar) + 1, length(p_bar) + 2))
sol <- csdp(C, A, b, K, control = csdp.control())
t(p_bar.mat) %*% sol$y
}
)
Error in library(plotly): there is no package called ‘plotly’ Traceback: 1. library(plotly) 2. stop(txt, domain = NA)
df <- data.frame(id = 1:length(optimal_res), return = optimal_res)
plot_ly(df, x = ~id, y = ~return, type="scatter", mode = "lines", height=500)
real.rets <- read.csv(
"~/Dropbox/Documents/Studies/MA (NRU HSE)/Year 1/Optimization Methods/Project/reveal_Opt/futures_returns.csv",
sep = ";", check.names = F, stringsAsFactors = F
)
real.rets <- real.rets[real.rets$Date > "2005-01-01", ]
real.rets <- real.rets[, -1]
real.rets <- real.rets
p_bar <- colMeans(real.rets)
real.rets <- real.rets[, which(is.finite(p_bar))]
p_bar <- colMeans(real.rets)
Sigma <- cov(real.rets)
Sigma_sqrt <- sqrtm(Sigma)
p_bar <- p_bar * 100
Sigma_sqrt <- Sigma_sqrt * 100
optimal_res <- sapply(
seq(0.0001, 0.1, 0.001),
function(eta) {
p_bar.mat <- matrix(p_bar)
C <- list(
matrix(0, length(p_bar) + 1, length(p_bar)+ 1),
c(rep(0, length(p_bar)), 1, -1)
)
A <- lapply(
1:length(p_bar),
function(i) {
x <- rep(0, length(p_bar))
x[i] <- 1
x <- matrix(x)
this.mat <- cbind(
rbind(
as.numeric(1 / qnorm(1 - eta) * t(p_bar.mat) %*% x) * diag(length(p_bar)),
t(Sigma_sqrt %*% x)
),
rbind(
Sigma_sqrt %*% x,
1 / qnorm(1 - eta) * t(p_bar.mat) %*% x
)
)
this.vec <- rep(0, length(p_bar) + 2)
this.vec[i] <- 1
this.vec[length(p_bar) + 1:2] <- c(1, -1)
list(
this.mat,
this.vec
)
}
)
b <- (-1) * p_bar
K <- list(type = c("s", "l"), size = c(length(p_bar) + 1, length(p_bar) + 2))
sol <- csdp(C, A, b, K, control = csdp.control())
t(p_bar.mat) %*% sol$y
}
)
df <- data.frame(id = 1:length(optimal_res), return = optimal_res)
plot_ly(df, x = ~id, y = ~return, type="scatter", mode = "lines", height=500)